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ABSTRACT 



Context. Thermal instability is one of the dynamical agents for turbulence in the diffuse interstellar medium, where both turbulence 
and thermal instability interact in a highly non-linear manner. 

Aims. We study basic properties of turbulence in thermally bistable gas for variable simulation parameters. The resulting cold gas 
fractions can be applied as parameterisation in simulations on galactic scales. 

Methods. Turbulent flow is induced on large scales by means of compressive stochastic forcing in a periodic box. The compressible 
Euler equations with constant UV heating and a parameterised cooling function are solved on uniform grids. We investigate several 
values of the mean density of the gas and different magnitudes of the forcing. For comparison with other numerical studies, solenoidal 
forcing is applied as well. 

Results. After a transient phase, we observe that a state of statistically stationary turbulence is approached. Compressive forcing 
generally produces a two-phase medium, with a decreasing fraction of cold gas for increasing forcing strength. This behaviour can 
be explained on the basis of turbulent mixing. We also find power-law tails of probability density functions of the gas density in 
high-resolution runs. Solenoidal forcing, on the other hand, appears to prevent the evolution into a two-phase-medium for certain 
parameter regions. 

Conclusions. The dynamics of thermally bistable turbulence show a substantial sensitivity to the initial state and the forcing properties. 
Key words, hydrodynamics - instabilities - turbulence - methods: numerical - ISM:kinematics and dynamics 



1. Introduction 

^ Observations of the interstellar medium (ISM) show large vari- 
. ations in the thermodynamical properti es of the gas, which le d 
to the famous three-phase-model of McKee & Ostriker (1977). 
, In this article, we focus on the warm neutral medium (WNM) 

■ and the cold neutral medium (CNM), wh ose occurrence was 
first explained consistently by iFieldl (1 19651) with the help of the 
thermal instability (TI). The classical conception of a more or 

■ less static multiphase medium with different phases coexisting 
] in pressure equilibrium (see, for instance, Wolfire et al.l 1 19951) 

, has nowad ays been replaced by a picture of a highly turbulent 

■ ISM (e. g. iBalle steros-P aredes et all 120071: iMcKee & Qstriked 
" l2007l) . Several observational results point to turbulence in the 

ISM, for example, the broadening of spectral lines, chemical 
mix ing of the ISM, or the star f ormation rate (for an over view, 
see iMac Low & KlessenI l2004l: lElmegreen & Scald l2004l) . As 
turbulence d ecays roughl y on on e dynamical timescale (e. g. 
iMac Low e t al. 1998; Stone et al.1'1998), energy has to be in- 
jected into the ISM to maintain turbulent motion. Several differ- 
ent sources were proposed to drive turbulence in the ISM, e. g., 
blast waves from supernova explosions, the differential rotation 
of galaxies, but also thermal or gravitational instabilities. 

While the properties of forced turbulence in the molecular 
and nearly isothermal phase have been studied extensively 
fe^ g. Krits uk et al.1 120071: ISchmidt et"aD l2009l: Federrath et al.' 
[2009a), there is no corresponding analysis for the diffusive 
ISM so far. Previous work on turbulence i n thermally bistabl e 
gas encompass decaying turbulence (.Kritsuk & Normanll2002ll) , 



simulations on kpc scales to study the effec t of stellar feedback 
(iJoung & Mac Lowl2QQ6l:lJoung et al.l20(39b . and c olliding flows 



of ato mic gas subject to the thermal instability (Heitsch et al 
120061 : [Vazquez-Semadeni et al. 2007; Hennebelle & Audit 
2007; BanerieeelaD [20091: lAudit & Hennebelld 120 icf 
FolinietaL .2009) . The last scenario has attracted a lot of 
attention as a possible model for the formation of molecular 
clouds from the diffuse ISM. Furthermore, simulations of 
turbulence produced by purely solen oidal forcing in thermally 



bistab le gas were performed by IVazquez-Semadeni et~aD 



(I2003h. iG azol et al.l ^ (i20051 iKissmann et al.l (12008"). and 
iGazol et al. ( 200^7 In these simulations, large amounts of 
thermally unstable gas (up to 50%) were found, which is also 
indicated by observational s tudies ( e. g. Dickey et al.l 1 1977 



Spitzer & Fitzpatrickl [T995I : iHeilesI 12001 



20031: iKanekar et al.ll2003h . 



Heiles & Trolandl 



As was demonstrated by ISchmidt et al.l (l2009l) for the 
isothermal case, compressive stochastic forcing produces a flow 
structure on large scales that resembles colliding flows, because 
converging shocks are formed between rarefied gas regions. In 
this work, we study the properties of compressively driven tur- 
bulence in thermally bistable gas, depending on the magnitude 
of the force and the initial state of the gas, which is given 
by the equilibrium temperature for the mean gas density. In 
this regard, our approach is similar to the work of iGazol et al.l 
(2005). The mean density is constant because of the periodic 
bounda ry conditions. The c o oling func tion follows the defini- 
tion of I Audit & Hennebelld (l2005h . As iFederrath et aP (l2008h 
reported a substantial sensitivity of the density field in isothermal 
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supers onic turbulence on the forcing (see also iFederrath et alJ 
l2009allbl ). we applied solenoidal forcing as well. Because of the 
three-dimensional parameter space (mean density, magnitude, 
and ratio of solenoidal to compressive modes of the forcing), 
the aim of the present study is to investigate the general be- 
haviour of forced thermally bistable turbulence by means of a 
series of moderate-resolution simulations, using uniform grids 
with 256^ cells. To avoid a further increase in dimensionality of 
the parameter space, we omit magnetic fields, self-gravity, and 
thermal conduction. Nevertheless, our results serve as a useful 
guide estimating the influen ce of unresolved thermal processes 
in galactic disk simulations ("Tasker & Bryan|[2QQ6l: iDobbs et^ 
i2008; Agertz et al. 2009). In such simulations, the cutoff length 
at the highest level of refinement is typically of the same order of 
magnitude as the forcing length scale in the periodic-box simu- 
lations. In the various prescriptions of the local star formation ef- 
ficiency, it is neglected that the fraction of cold gas that can turn 
into stars via gravitational collapse is less than unity. We there- 
fore investigate the influence of the forcing and the mean density 
on this cold gas fraction in detail. Future work will concentrate 
on particular parameter sets to obtain the turbulent scaling and 
fragmentation properties from high-resolution simulations and 
to study the influence of additional physics. 

The outline of this article is as follows. Section [2] briefly de- 
scribes the numerical techniques. The parameters, for which we 
run simulations, are specified in Sect. O The phenomenology of 
turbulence in thermally bistable gas is described for a represen- 
tative parameter set in Sect. IH For this parameter set, the grid 
resolution was varied from 128^ to 512^ to check for the resolu- 
tion dependency of the results. Next, we analyse the variation of 
statistical properties for diff'erent simulation parameters. In Sect. 
[3 the results are discussed in the context of other numerical stud- 
ies and observations, followed by a summary of the main results 
and an outlook in Sect. [6j 



2. Numerical techniques 

The simulation present ed in this article are performed with the 
open source code Enzo ( Q'Sh ea et al.l l2005). We solve the com- 
pressible Euler equations using the pi ecewise-parabolic method 
(PPM) of lColella & WoodwardI (11984 . Besides an external driv- 
ing force /, which is responsible for generating turbulent mo- 
tions, we include diff'erent heating and cooling terms, combined 
into a net cooling rate per unit mass X. The resulting equations 
for the mass density p, the velocity m, and the specific total en- 
ergy e of the fluid can be written as 



D 

— p = -pV • u , 



D 



D 



p—e + V 'Pu = pf ' u - p£(p,T) , 
where the total time derivative is defined by 



D 
Dt 



dt 



-hM • V. 



The specific total energy e is given by 



1 2 P 

2 (r-i)p 



(1) 
(2) 
(3) 

(4) 

(5) 



where y = 5/3 is the adiabatic exponent of monatomic gas, and 
the gas pressure P is related to the mass density and the temper- 
ature T via the ideal gas law: 



JUflH 



(6) 



The constants ^g, yu, and mn denote the Boltzmann constant, the 
mean molecular weight, and the mass of a hydrogen atom, re- 
spectively. 

The cooling function p£ was defined by lAudit & Hennebelld 

(12005 ). It includes the fine-structure cooling of Cii and Oi, as 
well as the cooling by H (Lya-line) and by electron recombi- 
nation onto charged grains. The only heating process consid- 
ered is the photoelectric eff'ect on small grains and polycyclic 
aromatic hydrocarbons (PAH) caused by far-ultraviolet galactic 
background radi ation. For deta iled information about the diff'er- 
ent processes seelWolfire et al.l (fl995. 2003); Spitzer (1978) and 
Bakes & TielensI ((1994). The cooling is implemented in the form 
of an explicit scheme into Enzo ( Nikl aus et al...2009) . After each 
time step, the state variables are updated in several subcycles by 
adding the resulting total energy increment. As the considered 
cooling and heating processes are only well defined in the diff'use 
ISM, the calculation of £(p,T) is constrained to temperatures T 
G [10 K, 10000 K]. For gas with temperatures above 10000 K, 
we use X(p,10000 K), which leads to an underestimate of the 
emitted energy, hence to warmer gas than expected. In strongly 
driven simulations we apply a temperature cutoff' at 30 000 K, 
assuming that the gas would cool down very fast beyond this 
limit. Since the speed of sound increases with the temperature, 
the cutoff' allows us to avoid very small time steps arising from 
unphysical values of the temperature outside the range of the 
cooling function. The concerned simulation runs are mentioned 
in Sect. II 

The specific force / in Eqs. ^ and ([3]), which generates 
large-scale motions of the gas, is composed in Fourier space. 
Each mode of the force field is given by a stoc hastic process that 
is based on t he Qrnstein-U hlenbeck process (lEswaran & Popd 
1988; Schmi dt et al.ll2006h . For the details of this method, see 
Sch midt et all (120091) . The forcing varies on the autocorrelation 
timescale T, which is given by the ratio of the integral length L 
to the characteristic velocity V of the large-scale motions. We set 
the integral length to half of the box size. Thus, the wavelengths 
of the forcing modes range from L/2 to 2L. The free parameter 
V specifies the magnitude of the forcing, i. e., the specific force 
is ~ V^/L. By applying a Helmholtz decomposition in Fourier 
space, the ratio of solenoidal (divergence-free) to compressive 
(rotation-free) modes can be arbitrarily adjusted. We use the de- 
composition parameter f such that ^ = I corresponds to purely 
solenoidal forcing, and the force becomes increasingly compres- 
sive as f decreases from 1 to 0. 



3. Simulation parameters 

Since we are interested in the local structure of the diff'use ISM, 
we treat the simulation domain as a part of a larger region with - 
in the statistical sense - homogeneous properties. Numerically, 
this is done by using a cubic simulation domain with periodic 
boundaries and a linear size of X = 40 pc. At the simulation start, 
the box is filled with atomic hydrogen of uniform density po and 
temperature Tq. For the mean particle densities no (which re- 
mains constant because of the periodic boundaries), we chose 
1.0 cm"^, 1.8 cm"^ and 3.0 cm"^, which are comparable to the 
mean density of the diff'use ISM. In the initial phase, the nearly 
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Fig. 1: Pressure equilibrium curve for the used cooli ng and 
heatin g processes. Via the isobaric instability criterion (iHuntej 
11970") the regions of thermal instability (red) and thermal stabil- 
ity (green) are determined. In addition, the isothermal lines at 
T - 200 K and T - 5000 K are shown, which define our region 
of thermally unstable gas. The dots represent the starting points 
of the simulations, which all lay in the thermally unstable region. 



homogeneous gas will adjust very fast to the equilibrium tem- 
perature that is given by the condition X = 0, regardless of the 
choice of Tq. For this reason, we set To equal to the equilibrium 
temperature corresponding to ^o- 

The pressure equilibrium curve Pcqin) following from £ = 
and the ideal gas equation (see Sect. [2]), is plo t ted in Fig. 
[TJ Via the isobaric instability criterion of Hunted (Il970l) , one 
can determine the regions of the thermal instability (TI) in the 
phase diagram that is shown in Fig. [TJ We define three diff'er- 
ent phases, a cold phase with T < 200 K, an unstable phase 
with 200 K < r < 5000 K, and a warm phase with T > 5 000 K 
This classification has also be en used in o ther works ( HeilesI 
120011 : [keiles & Troland 2003; lAudit & He nnebelle 2005^ and 
is warranted by the positions of the isothermal lines at 
T = 200 K, 5000 K, which - of course only very roughly - de- 
scribe the region of TI in the intermediate and high-density 
range. In the following, the terms cold, unstable, and warm phase 
always refer to the above definitions. The chosen initial con- 
ditions of our simulations are located in the thermally unsta- 
ble region (Fig.[T]). In consequence, one can expect that a two- 
phase-medium will evolve. As explained in Sect. O we start the 
simulations with the gas at rest, and turbulence is generated by 
stochastic forcing. The forcing strength is regulated by the char- 
acteristic velocity V. Normalising V by the initial soundspeed 
^^0 = (ykBT{) I iiniHy^ , we obtain a characteristic Mach number. 
Ma = y/co, which determines the magnitude of the stochastic 
force. We performed simulations with Ma between 0.2 and 3.0. 
However, Ma is not an indicator for the average Mach number 
of the flow, which depends on the phase separation of the gas. 

In one set of simulations, we applied purely compressive 
forcing (^ = 0) for difl'erent initial densities. In addition, simula- 
tions with ^0 = 1-0 cm"^ were also carried out for the solenoidal 
case (f = 1). The parameters of all simulations are listed in Table 
[T] Each simulation ran until an approximate statistically steady 
state was found for at least one dynamical timescale T. Although 
it is not obvious that a well-defined statistically steady state can 
be reached for non-isothermal turbulence, all simulations pre- 
sented here show a convergence of statistical quantities related 



Table 1: Performed simulations with the mean density hq, the 
starting temperature To, the Mach number Ma, f and the simu- 
lation time, normalised to the dynamical timescale T. 



Run 


riQ (cm"^) 


ro(K) 


Ma 




i{T) 


DlMaO.2 


1.0 


4908 


0.25 





6.0 


DlMaO.S 


1.0 


4908 


0.5 





7.5 


DlMal.O 


1.0 


4908 


1.0 





6.0 


DlMaL5 


1.0 


4908 


1.5 





5.0 


DlMa2.0 


1.0 


4908 


2.0 





5.0 


DlMaL0_128 


1.0 


4908 


1.0 





6.0 


DlMal.0_512 


1.0 


4908 


1.0 





6.0 


D1.8Ma0.2 


1.8 


1353 


0.2 





4.0 


DL8Ma0.5 


1.8 


1353 


0.5 





5.0 


Dl.SMal.O 


1.8 


1353 


1.0 





5.0 


DL8Ma2.0 


1.8 


1353 


2.0 





5.0 


D1.8Ma3.0 


1.8 


1353 


3.0 





4.0 


D3Ma0.5 


3.0 


493 


0.5 





3.0 


D3Mal.O 


3.0 


493 


1.0 





5.0 


D3Ma2.0 


3.0 


493 


2.0 





5.0 


D3Ma3.0 


3.0 


493 


3.0 





5.0 


DlMa0.5s 


1.0 


4908 


0.5 


1 


3.0 


DlMal.Os 


1.0 


4908 


1.0 


1 


5.0 


DlMa2.0s 


1.0 


4908 


2.0 


1 


5.0 



to the flow properties and the thermodynamic properties of the 
gas as time proceeds. Thus, the statistics described in the next 
sections are averaged over the last dynamical timescale. Since 
the computational cost of performing simulations of turbulence 
in non-isothermal gas is significantly higher than for isothermal 
gas, all runs were performed on a static grid of N = 256^ cells. 
The only exceptions are the runs DlMal.0_128 with a resolution 
of 128^ and DlMal.0_512 with 512^ to test the resolution de- 
pendency of the results. Because of the moderate resolution of 
the simulations, we only investigated global statistics and prob- 
ability density functions, but not two-point statistics or power 
spectra. 

4. Results 

4.1. Basic properties and resolution study 

A crucial question for the following analysis is to what extent 
statistical properties are aff'ected by the numerical resolution. To 
investigate the dependence on resolution, we consider the simu- 
lations DlMal.0_128, DlMal.O and DlMal.0_512 with com- 
pressive forcing, a mean number density no = 1.0 cm"^, and 
Ma = 1.0. The resolution varies from 128^, over 256^ to 512^. 
Figure [2] shows the time evolution of the temperature field and 
the Mach number for these three simulations. For them, we see 
a similar evolution in time both for the temperature field and 
the velocity field (represented here by the Mach number). The 
onset of the TI becomes apparent by a rapidly spreading tem- 
perature distribution within the first dynamical timescale, where 
the temperatures finally range from roughly 10 K to several 
10 000 K. For the lower resolutions, we did not apply the cutoff* 
at 30 000 K. As one can see, the temperature rises locally to as 
much as 10^ K. For the simulation with 512^ grid cells, these ex- 
trema entailed problems with the time stepping. Consequently, 
we introduced the high-temperature cutoff* for t > 3 T. After 
a few dynamical timescales a statistically steady state is ap- 
proached, in which the mass- and volume- weighted mean values 
are nearly constant. As can be seen in Fig.|2l the mass-weighted 
mean of the temperature, T'mw, is much lower than the volume- 
weighted temperature T'vw, because the cold gas is associated 
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(e) DlMaL0_512 (f) DlMaL0_512 

Fig. 2: Time evolution of the temperature field (left panel) and the Mach number (right panel) for the simulation runs DlMal.0_128 
(upper panel), DlMal.O (mid panel), and DlMal.0_512 (lower panel). Plotted are the volume- weighted mean values, the mass- 
weighted mean values, and the maximum and minimum values. 



with high-mass densities. The low temperature of the cold, dense 
gas results in a supersonic mass-weighted mean Mach number, 
Almw, while Alvw is transonic. 

Averaging the diff'erent mean values over the last dynami- 
cal timescale, ()^, we only find small diff'erences between the 
time-averaged mean values for the diff'erent resolutions. For 
the mass-weighted temperature mean values, which correspond 
to the mean thermal gas energy, we obtain (7~mw)? = 2229 K, 
2333 K, and 2326 K for simulation DlMal.0_128, DlMal.O, 



and DlMal.0_512, respectively. Thus, the relative deviations are 
less than 5%. For the volume- weighted values <7~vw)?, the devi- 
ations are even less than 0.6%. For (Alvw)? and (Atvw)?, on the 
other hand, the diff'erences are below 3.7% and 20%, respec- 
tively. The larger deviations of the time-averaged mean Mach 
numbers point toward a stronger sensitivity of the turbulent ve- 
locity field on the numerical resolution, while the thermody- 
namic properties appear to be quite robust. 
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Fig. 3: Left: Temperature PDFs for the runs DlMal.0_128, DlMal.O and DlMal.0_5 12, averaged over the last dynamical timescale. 
The maxima around 55 K and 6500 K indicate the CNM and the WNM. Right: Density PDFs for the same runs as in the left panel. 
For comparison a power law with an exponent of -3.2 is shown. 



Further indicators for the resolution dependency are proba- 
bility density functions (PDFs). The PDFs of the number density 
and the temperature, averaged over the last dynamical timescale 
in the simulations DlMal.0_128, DlMal.O, and DlMal.0_512, 
are plotted in Fig.O The two maxima of the temperature PDFs 
show that a cold and a warm gas phase are formed, with some 
gas in the unstable phase in between. The positions of the two 
maxima around 55 K and 6500 K are almost independent of 
the numerical resolution. These temperatures can be associated 
with the CNM and the WNM. Near the extrema of the temper- 
ature, there are slight deviations between the PDFs for differ- 
ent resolutions. The small excess of low-temperature gas in run 
DlMal.0_512 corresponds to a slightly higher amount of very 
dense gas. The tiny fraction of gas hotter than 10 000 K supports 
our assumption that the cutoff at 30000 K in run DlMal.0_512 
has no significant impact on the dynamics of the whole system. 
In the temperature interval from roughly 100 K up to 1000 K, 
we see a decrease in the amount of gas as the resolution in- 
creases. Consequently, the phase separation becomes more pro- 
nounced at higher resolution. The density PDFs shown in the 
right panel of Fig. [3] reveal a small excess of high-density gas 
for simulation DlMal.0_512, which corresponds to the afore- 
mentioned increase in cold gas. For the highest resolution, the 
density PDFs also show a decrease in the amount of gas in 
the intermediate density range, which corresponds to thermally 
unstable gas. A notable feature in the high-density range is 
the power-law shape, which extends over more than one or- 
der of magnitude. We estimate the power-law exponent of run 
DlMal.0_512 to roughly -3.2, whereas no power-law exponents 
are calculated for the lower resolution runs. Nevertheless, the 
fit for run DlMal.0_512 seems to apply quite well even for run 
DlMal.0_128 and DlMal.O, as shown in the right panel of Fig. 
[21 The power law at these high densities is caused by an effective 
equation of state (EOS) in this range, which shows a polytropic 
index yeff < 1 • 

Although the statistics we investigated indicate that the simu- 
lations DlMal .0_128 and DlMal .0 are not fully converged with 
respect to the density and temperature distributions, we conclude 
that the approximate behaviour of thermally bistable turbulence 
can be inferred from simulations with 256^ grid cells. This is the 
aim of the parameter study that is presented in the following. 

To illustrate the properties of thermally bistable gas. Fig. |4] 
shows contour plots of the particle density, the temperature, the 
vorticity norm and the Mach number in 2D slices for simula- 



tion DlMal.0_512 in the statistically stationary regime. When 
comparing the slices in the top panel of Fig.lH one can see that 
regions of high gas density correspond to low temperatures and 
vice versa. The most prominent feature that can be seen in the 2D 
slices is the big clump of cold, dense gas. While the boundaries 
of this clump are rather sharp, there are also regions in which 
the cold and warm gas phases entrain each other and form intri- 
cate structures. The rich small-scale structure becomes apparent 
in projections of the mass density and the temperature, which 
are shown in Fig. \5\ Even in these projections, relatively large 
clumps of cold gas appear prominently. Remarkably, the vor- 
ticity is largely reduced in these clumps in comparison to the 
surrounding medium, which is filled by more or less homoge- 
neous turbulence (see left bottom panel of Fig.|4]). Nevertheless, 
the two-dimensional PDF of the vorticity magnitude versus the 
density, which is shown in Fig. [6l demonstrates that high vor- 
ticity is found for the whole range of gas densities. Also the 
density-dependend mean value shows no clear trend but is more 
or less constant over the whole range of densities. This seems to 
contradict the observation made in the vorticity slice in Fig. |4l 
A possible explanation could be that the vorticity remains low 
only in larger clumps of overdense gas, while in smaller clumps 
or in thin filamentary, but nevertheless overdense structures, the 
vorticity is about as high as in the surrounding low-density gas. 
Presumably, turbulent eddies that are produced by the TI in the 
unstable phase rapidly decay inside compact cold gas regions, 
because of the short sound crossing time. The compressive forc- 
ing, on the other hand, does not directly generate vorticity and it 
fails to induce shocks inside the clumps of extremely dense gas. 
It is possible, however, that these clumps would fragment if the 
numerical resolution were increased. The Mach number, on the 
other hand, indicates supersonic motions in the cold, dense gas 
and subsonic flow in the WNM (see right bottom panel of Fig. 
[4]). This leads to the effect that the mass weighted mean value 
Almw is higher than the volume weighted mean value Alvw as 
shown in the right panel of Fig.O 

Slices of the gas density and the vorticity after 1 and 2 in- 
tegral timescales, which are shown in Fig. [71 demonstrate that 
turbulence is very inhomogeneous at these stages of the flow 
evolution. The patches of high vorticity that can be seen are as- 
sociated with density fluctuations and shocks in the unstable or 
warm gas phases. Thus, it appears that the initial production of 
vorticity mainly results from the forcing. As the TI becomes ac- 
tive and temperature fluctuations grow, this production mecha- 
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Fig. 4: Temperature (top left), density (top right), vorticity field (bottom left), and Mach number (bottom right) in a slice through 
the simulation domain of run DlMal.0_512 at t = 6.0 T. 



nism is accompanied by baroclinic vorticity generation (see, for 
instance Schmidt 2007). In the late, nearly stationary phase, the 
high vorticity is found for the whole range of gas densities (see 
Fig.O. 



4.2. Influence of the forcing strength 

So far, we have considered a characteristic Mach number 
Ma = 1, which means that the velocity of the flow produced 
by the forcing is close to the speed of sound for the initial 
gas temperature. Next we analyse the influence of varying the 
strength of compressive forcing for no = 1.0 cm^ and N = 256^. 
For the runs DlMal.5 and DlMa2.0, the temperature cutofl' at 
30 000 K has been applied to avoid artificially high tempera- 
tures. In Fig. [8l the temperature and density PDFs are shown 
for the characteristic Mach numbers in the range from 0.2 to 
2.0. In each case, the PDFs are averaged over the last dynam- 
ical timescale of the simulation. Except for the temperature 
PDF of run DlMa2.0 all other PDFs in the left panel of Fig. 
[8] show two maxima located in the temperature range of the 
CNM and the WNM. The positions of the peaks do not vary 
significantly with the forcing strength in contrast to the peak 
width. The clearly visible bimodal temperature distribution for 



the weakly driven simulations DlMaO.2 and DlMaO.5 gradu- 
ally vanishes with increasing forcing strength, while the amount 
of gas in the unstable range increases. In run DlMa2.0, no bi- 
modal temperature distribution can be discerned. The density 
PDFs of the weakly driven simulations with Ma < 1 in the right 
panel of Fig. [8] also show a bimodal shape, where the strong 
peak in the low-density range (~ 0.2 cm"^) and the plateau 
around 10 cm"^ correspond to high- and low-temperature peaks, 
respectively. With increasing forcing strength, the maximum of 
the density PDF broadens substantially. Thus, for high charac- 
teristic Mach numbers, the density distribution appears to be 
increasingly dominated by the nearly adiabatic expansion and 
compression waves produced by the forcing. Moreover, the high- 
density tails become flatter; i. e., stronger forcing tends to induce 
higher densities. However, even for DlMa2.0, the peak densities 
up to several 1000 cm"^ are much higher than what is observed 
in isothe rmal simulations wit h comparable characteristic Mach 
number (ISchmidt et al.l l2009b . Consequently , the high-density 
tails naainly results from the TI (see also Kritsuk & NormanI 
2002^; iLi etal.1 [2003b . Accordingly, all density PDFs exhibit a 
power-law like shape in the high-density range, although the nu- 
merical resolution is too low to draw definite conclusions. 
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Fig. 6: Two-dimensional PDF of the vorticity vs. the density for 
run DlMal.0_512 at ^ = 6.0 T with logarithmic scaled contour 
lines. The thick black line shows the density-dependent mean 
vorticity. 



The pressure-density phase diagrams for DlMaO.2 and 
DlMal.O, which are shown in Fig.O demonstrate how the ther- 
modynamic state of the gas is affected by turbulence. In the 
case of turbulence with low intensity (left panel), the distribu- 
tion shows t ypical properties o f a static two-phase-medium as 
described bv lField et al.] ([T969h . Here the gas exists in a rough 
pressure equilibrium with deviations from the mean pressure less 
than half an order of magnitude. Only the high-density gas is out 
of pressure equilibrium, because of the very small cooling times 
corresponding to these densities. For stronger turbulent motions, 
as in run Dl Mai .0 (right panel of Fig. O, the gas deviates largely 
from pressure equilibrium. Thus, for gas in a strongly turbulent 
state, the classical idea of two distinct phases, coexisting in a 
rough pressure equilibrium with only small amounts of gas in 
the unstable range ceases to apply. 

To quantify the increasing amount of thermally unstable gas 
that is observed for strong forcing amplitudes, we determined the 
mass fractions of the warm, unstable and cold phases as defined 
in Sect. O The results are shown in Fig. [TOl where the diff'er- 
ent mass fractions are plotted against the total specific kinetic 



energy for each simulation. In accordance with the PDFs, the 
mass fractions, as well as the specific kinetic energies, are aver- 
aged over the last dynamical timescale. Clearly, the turbulence 
energy resulting from the forcing has a strong influence on the 
distribution of the gas into the diff'erent phases. The mass frac- 
tion of the unstable gas increases at the cost of the mass fraction 
of the cold gas, whereas the mass fraction of the warm phase 
is only weakly afl'ected. That the redistribution of the gas with 
increasing forcing strength mainly occurs between the cold and 
the unstable phase and, to a lesser extent, between the warm and 
the unstable phase can be understood as the consequence of the 
low compressibility of the hot gas, in which only sub- or tran- 
sonic motions are induced, regardless of the characteristic Mach 
number. The dynamics in the unstable phase and the cold phase, 
on the other hand, are highly sensitive to the forcing strength. 
Although compressive forcing produces the seeds for the TI, we 
note that high turbulence intensity inhibits the gas from settling 
into the cold phase. However, it appears that the TI overcomes 
this efl'ect once a critical amount of cold gas is accumulated. 
This is suggested by the occurrence of compact regions of cold 
gas with weakly turbulent interiors (see Fig.©. 

4.3. Tracer particle statistics 

To examine the cause of the higher mass fraction of gas 
in the thermally unstable range for stronger forcing, we in- 
cluded 150000 Lagrangian tracer particles in the simula- 
tions DlMaO.2 and DlMa l.O. This approach is motivated by 
lAudit & HennebeUd (l2005h . who also found an increase in the 
thermally unstable gas with increasing turbulence intensity at the 
cost of the cold gas mass. The authors explain this efl'ect on the 
basis of a simple semi-analytical model. This model implies that 
the time interval that is spent by a fluid element in the unsta- 
ble phase before returning to the stable phase becomes longer as 
the turbulence intensity increases. On the other hand, the mass 
of gas in the unstable phase is also controlled by the frequency 
of strong perturbations that bring the gas out of the equilibrium 
state. Thus, the larger amount of gas in the unstable phase can 
also be interpreted as a consequence of enhanced turbulent mix- 
ing. With the help of the tracer particles, we are able to test the 
predictions of this model by explicitly following the time evo- 
lution of a fluid element moving with the flow. In particular, we 
measure the time-interval ^up,/, defined as the time, a tracer parti- 
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cle is located in the unstable range before returning to the stable 
state of the CNM or the WNM. We find that the mean duration 
^up of the stay of a tracer particle in the unstable phase is re- 
duced by a factor of 2.3 if the characteristic Mach number of 
the forcing increases from 0.2 to 1 .0. However, after analysing 
how often the particle is perturbed out of the equilibrium state 
into the unstable phase, we find that the frequency of the pertur- 
bations, , rises by a factor of ~ 4.5 from run DlMaO.2 to 
DlMal.O. This implies that the thermally unstable gas fraction 
roughly increases from run DlMaO.2 to DlMal.O by a factor of 
4.5 • 23 ~ 2.0. 

4.4. Variation in the mean density 

It is also of interest to examine how the gas distribution changes 
if the mean density of the gas is altered. Interestingly, it turns out 
that even a small increase of the mean density has a noticeable 
eff'ect on the distribution of the gas among the diff'erent phases. 
In particular, we performed several simulations with compres- 
sive forcing for mean densities of 1.8 cm"^ and 3.0 cm"^ (see 
Table [B. Figure [T] shows that the initial state of the gas in these 
simulations is located closer to the local minimum of the equi- 
librium curve. For = 3.0 cm"^, the initial state is located near 
the low-temperature boundary of the central unstable region. The 



dependence of the mass fractions of the diff'erent gas phases on 
the forcing strength is shown in Fig. [TTJ While the qualitative 
behaviour is similar to what is observed for no = 1.0 cm"^ (see 
Fig. [T0|), the fraction of cold gas is clearly greater for a higher 
mean density. In the case = 3.0 cm"^, about 90 % of the mass 
settles into the cold phase, more or less independent of the tur- 
bulence energy. This result is reasonable because, as the initial 
state is chosen closer to the cold phase, an increasing amount 
of the gas will cool and settle into the cold phase due to the 
external compressions. At densities no > 10 cm^, most of the 
gas resides in the nearly isothermal right branch of the equilib- 
rium curve, and only strong rarefactions can push the gas into the 
warm phase. In this regime, the dynamical properties are com- 
parable to those of isothermal turbulence. 

4.5. Solenoidal forcing 

ISchmidt et all (l2009l) and iFederrath et al.l (l2009ah show in their 
analysis of supersonic isothermal turbulence that there are sig- 
nificant diff'erences between compressive and solenoidal forc- 
ing. In particular, the density PDFs are m uch broader for com- 
pressively driven turbulence (iFederrath e t al. 2008, 2009a). For 
non-isothermal gas, the diff'erences should be even more pro- 
nounced because the TI is aff'ected by the diff'erent density dis- 
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tributions, which are generated by distinct forms of the forcing. 
The PDFs of the number density and the temperature result- 
ing from three simulations with solenoidal forcing and a mean 
density of 1.0 cm"^ (listed in Table [B are plotted in Fig. [121 
where the temperature cutoff at 30 000 K was applied for the 
runs DlMal.Os and DlMa2.0s. As one can see, there is no phase 
separation at all. Both the density and temperature PDFs show a 
single peak that roughly corresponds to the initial state. With 
increasing forcing strength, the distributions become broader, 
and the gas tends to become hotter. The higher temperatures 
that are found for Ma = 2.0 in comparison to the other runs 
are caused by the stronger dissipative heating that balances the 
higher rate of energy injection by the forcing. The amount of 
cold gas with temperatures below 200 K is zero. Thus, it ap- 
pears that the TI is largely suppressed by solenoidal forcing, 
although a significant fraction of the gas is in the unstable re- 
gio n (see Fig.[T])- The density PDFs are similar to those reported 
by iPassot & Vazquez-Semadenil (1 19981) for simulations of tur- 
bulence in a polytropic gas with y > 1, particularly, with re- 
gard to the low-density tails. At least for the runs with Ma < 1 .0, 
there appears to be a power-law range for densities lower than 
1.0 cm-^ 
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Fig. 10: Time-averaged mass fractions of the warm (red points), 
unstable (green points), and cold phase (blue points) for simula- 
tions with ^0 = 1-0 cm"^ and N = 256^ plotted against the spe- 
cific kinetic energy. An increasing kinetic energy corresponds to 
an increase in the forcing Mach number (compare Table [TJ. 



5. Discussion 

Most simulations considered in this article were run with a res- 
olution of N = 256^ grid cells. For the chosen box size of 
40 pc, the cutoff' length in these simulations is A = 0.16 pc. 
Consequently, cold gas structures with a typical size in the range 
of a few parsecs are only marginally resolved and substantially 
smeared out by numerical diff'usion. For this reason, it cannot be 
expected that the properties of these structures can be described 
in a detailed quantitative manner. However, the resolution study 
we performed for a particular parameter set suggests that the tur- 
bulent phase separation of the gas can be approximated reason- 
ably well. In particular, we showed that trends can be seen if the 
forcing is varied, which allows for a qualitative discussion. This 
will help us selecting parameters for more elaborate numerical 
studies at higher resolution, using adaptive mesh refinement. In 
particular, an enhanced numerical resolution will be necessary 
to confirm the indications for power-law tails of the mass den- 
sity PDFs toward high densities and to determine the slope with 
higher precision. These power laws can be found in nearly all 
simulations of compressively driven turbul ence in the prese nt 
study. This clearly supports the res ults of IScalo et al.l ( 1998i) . 
iPassot & Vazquez-Semadenil (Il998h . and iLi et al.1 (120031) . who 
applied a simple polytropic EOS with y < 1 as a simple model 
for cooling. These authors obtained skewed lognormal density 
PDFs with power-law tails at high densities. In our simulations, 
the thermodynamic state of the gas at high densities roughly fol- 
lows the equilibrium curve, which corresponds to an eff'ective 
polytropic index y^s < 1. Power-law tails for non-isothermal 
turbulen ce followed also from the s imulations of decaying turbu- 
lence by lKritsuk & Norni aiJ ('2002') and from a recent colliding- 
ffow simulation of lAudiTfe HennebeUi (l2010h . in which a poly- 
trope with y = 0.7 was applied. 

The main result of our study is the increase in the mass of the 
thermally unstable gas for higher turbulence energy (Fig. [TOl) , 
which is roughly balanced by a decrease in the mass of cold 
gas. Consequently, stronger compressive forcing does not result 
in more cold, dense gas. The reason is found in the broadening 
of the pressure distribution with increasing forcing strength and, 
thus, kinetic energy of the turbulent ffow. While weak forcing 
produces the initial seeds for the TI by perturbing the density 



of the gas, the forcing drives the gas away from a state of ap- 
proximate pressure equilibrium if the turbulent velocity becomes 
comparable to the speed of sound. This eff'ect is shown clearly 
for runs DlMaO.2 and DlMal.O in Fig.[9l Only at high densities 
is the gas mainly concentrated around the equilibrium curve, in- 
dependent of forcing strength. This is a consequence of the very 
short cooling time in comparison to the dynamical timescale and 
the sound crossing time in this range. At low densities, on the 
other hand, the gas in run DlMal.O is nearly adiabatic because 
of the long cooling time in the WNM. The eff'ective polytropic 
index in this regime is found to be yeff = 1.52. As one can see in 
Fig. [8] (right panel), this entails a pronounced broadening of the 
density PDFs in the un der-dense range if Ma > 1 , which is a gain 
similar to the results of iPassot & Vazquez-Semadenil (Il998h for 
their simulations with y > I. 

In Sect. 14.31 we show with the help of tracer particles that 
the higher frequency of the perturbations of the particles out 
of the stable regime, , overcompensates the decrease in the 
mean duration, ^up , a tracer particle stays in the unstable phase 
if the forcing strength increases. This results in an increase in 
the unstable mass fraction roughly by a factor of 2.0 for sim- 
ulation DlMal.O compared to DlMaO.2, which approximates 
the ratio of 2.7 that follows from the mass fractions of unsta- 
ble gas in these simulations reasonably well. Our observation 
that ^up decreases with increasing forcing strength is contrary to 
the assumption of Audit & Henn ebelle (2005) that ^up increases 
with the intensity of turbulence. The discrepancy could appear 
because their model is 2-dimensional and applies to the unper- 
turbed transition of a ffuid element from the WNM to the CNM, 
which might be insufficient to describe 3 -dimensional turbulence 
with external forc ing. In order to rule out tha t the discrepancy 
with the model o f lAudit & Hennebelld (l2005l) arises because it 
only describes transitions from the WNM to the CNM, we also 
examined phase transitions between the WNM and the CNM 
with the tracer particles. We found that the mean transition time 
%NM^CNM decreases from simulations DlMaO.2 to DlMal .0 by 
a factor of ~ 2.5, which is similar to the behaviour of ^up- The 
averaged time intervals ^cnm and %nm, a ffuid particle stays in 
the CNM or in the WNM, become signiffcantly shorter as well. 
Summarising, the tracer particle statistics imply that phase tran- 
sitions between the different phases occur more frequently with 
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Fig. 9: Pressure-density phase diagrams of run DlMaO.2 (left) and DlMal.O (right) with a logarithmic scaling of the contour lines. 
The equilibrium between heating and cooling is indicated by the solid black line. 



increasing forcing strength, i. e., turbulent mixing is more ef- 
fective, which results in a higher mass fraction of the unstable 
phase. 

An important question is for which simulation parameters 
the mass fractions of the different phases that are plotted in Figs. 
[TOland[TT]are comparabl e to observational results for the atomic 
ISM (see, for instance, Heiles"20()l|; iHeiles & Trolandll2QQ3h . 
One difficulty in such a comparison is that observational val- 
ues usually do not refer to a particular mean density of the gas, 
for which measurements are taken. Moreover, the observations 
presented by these authors do not belong to certain diffuse H i- 
clouds but are absorption measurements on cold gas along the 
lines of sight against different background radio sources. Besides 
the mean Mach number in the cold phase, which the authors de- 
termined to be around 3, particular properties of these objects are 
not known. Nevertheless, we note that Heiles & Troland (20031) 
find that the average mass fraction of atomic gas in the ISM with 
temperatures in the range from 500 K to 5000 K is at least 29 
%. Together with the warm gas fraction, their measurements im- 
ply that the mass fraction of cold gas should be less than 40 %, 
which is clearly at the lower end of the range of values follow- 
ing from our numerical study. Thus, strong compressive forcing 
of turbulence in gas at a mean density of 1.0 cm"^ appears to 



be typical of the conditions in the diffusive ISM. However, the 
sources studied by Heiles & Troland (2003) show a large scat- 
ter around the average values, and there exist regions with mass 
fractions of cold gas up to nearly 90%, which would correspond 
to weak forcing and a higher mean density in our study. 

A surprising result is that we did not find any cold phase at 
all in the simulations with purely solenoidal forcing, regardless 
of the magnitude of the force. Although the turbulent velocity 
fluctuations in compressible gas entail density fluctuations and 
the initial state of the gas is in the unstable regime, it appears 
that the phase separation is completely suppressed. One possible 
explanation is that the TI would set in on length scales that are 
smaller than the numerically resolved scales in our simulat ions. 
The minimum length scale of the TI is the Field length (IFieldl 
J 965), which is defined by 
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With this, we get a Field length in the CNM that is at least 2 
or 3 orders of magnitude smaller than the spatial grid resolu- 
tion. Therefore, a cold phase might develop at higher resolu- 
tion. iHennebelle & Audit! (l2QQ7l) show in a 2-dimensional res- 
olution study that a spatial resolution below 2 • 10"^ pc is nec- 
essary to reach numerical convergence for a colliding flow of 
thermally bistable gas, which is unfeasible for 3 -dimensional 
simulations with the current computing power. However, the 
observed behaviour could also be physical and result from the 
apparent diff'erences between the turbulent velocity field pro- 
duced by s olenoidal forcing conip ared to compressively driven 
turbulence ( Schmidt et al.l l200 8b*). There are indications of an 
extended range of length scales, for which rotationally driven 
turbulence has not only lower density fluctuatio ns but also ex- 
hibits a more space-filling vorticity field ( Feder rath et al. I [20081 
I2009ah . As a consequence, the turbulent pressure could stabilise 
the gas eff'ectively against the TI. An analogous eff'ect has been 
observed for forced self-gravitating turbulence in a n isothermal 
gas (e. g. lKlessen et al ] l2000l:ISchmidt et al.ll2008ah . In addition, 
it is likely that the shearing eff'ect of the velocity fluctuations 
produced by solenoidal forcing disrupts incipient small conden- 
sati ons of dense gas. On th e other hand, Gaz ol et aP (l2005ll20Q9l) 
and iKissmann et al.l (12008.) found a cold phase for turbulence 
in non-isothermal gas that was produced by solenoidal forc- 
ing. These findings do not necessarily contradict our results, be- 
cau se [Ki ssmann et al. (2008) used a mean density of 2.0 cm"^, 
and lGazol et al.1 (l2005Ll2009h appHed a sHghtly diff'erent cooHng 
function. Thus, the initial state was located deeper inside the un- 
stable regime in both cases. For this reason, the TI can overcome 
the inhibiting eff'ect of turbulence more easily than in our simu- 
lations with ^0 = 1-0 cm"^. Indeed, we also found a cold phase 
developing in simulations with solenoidal forcing if the mean 
density was set to 3.0 cm"^. Taking these results together, it ap- 
pears that the evolution of thermally unstable gas is extremely 
sensitive to the initial conditions and the details of the cooling 
function if the forcing is solenoidal. 

The strong sensitivity on the mean density and on the forc- 
ing shown here could in principle be used to draw some conclu- 
sions about the forcing modes in observed diff'use H i-clouds. In 
practice, however, it is hard to state something about the under- 
lying forcing because of the difficulty determining the mean gas 
density observationally, and hence deciding whether the gas is 
thermally unstable or not. Apart from this difficulty, this work 
suggests that, for a marginally thermally unstable cloud (^o - 1 
cm"^) showing a bimodal temperature distribution, the forcing 
should mainly be compressive. On the other hand, for clouds in 
the highly unstable regime (^o ^ 3 cm"^), it is hardly possible 
to distinguish between solenoidal or compressive forcing as both 
forcing modes would produce a two-phase-medium with gas dis- 
tributions that are mostly determined by the thermal instability. 



6. Conclusion 

In this paper, we have performed a parameter study of forced 
turbulence in thermally bistable gas, with a cooling function that 
was defined by Audit & Hen nebelld (I2005D . Considering both 
compressive and solenoidal forcing, we find substantial varia- 
tions in the global statistics and probability density functions of 
the gas. Because of the limited resolution of the simulations, we 
did not compute two-point statistics or power spectra, nor did 
we analyse fragmentation properties. In the following, we sum- 
marise our results. 



1. Although we clearly see a phase separation for compres- 
sively driven turbulence already at moderate numerical res- 
olution, the PDFs of the mass density and the temperature 
suggest that a resolution higher than 512^ is required to ac- 
curately determine the tails of the distributions and the mass 
fractions of the different phases. For a forcing length scale 
around 20 pc, a spatial resolution significantly lower than 
0.1 pc is implied. 

2. A similarity of the density PDFs of all compressively 
driven simulations is the power-law like shape in the high- 
density range, which is caused by an effective polytropic 
index yeff < 1- Exploring the high-density tails in simula- 
tions with enhanced numerical resolution will bear conse- 
quences on theoretical predic tions for the core mass function 
(iHennebeUe & Chabrie3[2QQ 9). 

3 . For weak forcing, the warm and the unstable phase are nearly 
isobaric. With increasing forcing strength, we find larger de- 
viations from the pressure equilibrium and a greater amount 
of gas with temperatures in the intermediate, unstable range, 
which is roughly balanced by less gas in the cold phase. 
The Lagrangian statistics inferred from tracer particle sug- 
gest that this behaviour is caused by stronger turbulent mix- 
ing. 

4. The distribution of the gas among the diff'erent phases also 
depends on the mean density. For mean number densities 
higher than 1.0 cm"^, most of the gas tends to settle into 
the cold phase. 

5. If the forcing is purely solenoidal, we do not find bimodal 
density and temperature distributions for a mean density of 
1.0 cm"^. In this case, all the gas is in the thermally unstable 
and warm phase. 



Qualitatively, our results agree with oth er numerical studies of 
turbulence in non- i sothermal gas (e. g . Kritsu k & Normanll2002l: 
iGazol et al.1 [20051: lAudit & Hennebelle 2005). The large-scale 
flow structure produced by compressive stochastic forcing re- 
sembles the colliding flow scenario, since transient frontal col- 
lisions of shocks occur randomly ( Schmidt et al. 2009). Given 
the sensitivity of the phase separation on the mean density, a 
crucial difference arises from the accumulation of mass in the 
colliding flow scenario, because of the inflow boundary condi- 
tions. In a periodic box, on the other hand, the mass is con- 
served. Consequently, in colliding flow simulations, there is 
an implicit drift through the parameter space in the course of 
time. In this regard, we note that .Audit & Hennebelld (1201 Ol) 
do not find clear indications of a power-law tail of the den- 
sity PDF in their simulations with cooling. If self-gravity is in- 
cluded, a caveat of periodic-box simulations is that the amount 
of mass that is confined in gravitationally collapsing objects 
cannot indefinitely increase at a constant rate. A net inflow of 
gas from the boundaries, on the other hand, constantly reju- 
venates the supply of gas that can be bound by gravitational 
collapse (e. g. Yazquez-Semadeni et aDl2007l: [Hennebelle et al.l 



|2008[ iBaneriee et al.l 12009^ although a statistically stationary 
state cannot be reached either because of the reasons men- 
tioned above. Whether the physical conditions in the interstellar 
medium on a length scale of a few 10 pc are approximated more 
closely in simulations of colliding flows without external forc- 
ing compared to forced turbulence simulations at a fixed mean 
density remains a matter of debate. Regardless of this, the re- 
sults presented here show a substantial sensitivity on the initial 
state of the gas and the forcing properties. This illustrates that 
a restriction to highly specialised conditions possibly leads to 



12 



1 

0.8 
0.6 
0.4 
0.2 




D. Seifried et aL: Forced turbulence in thermally bistable gas: A parameter study 

1 



cold phase 
unstable phase 
warm phase 



10" 



10 



10 



2, 2. 



10^ 



10 



12 



0.8 
0.6 
0.4 
0.2 




cold phase 
unstable phase 
warm phase 



10" 



10 



10 



ej^j^/(cm /s ) 



2, 2, 



10 



11 



10 



12 



ej^j^/(cm /s ) 
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misleading conclusions, as even minor changes in the studied 
conditions can cause significantly diff'erent results. 

In large-scale s i mulat i ons of the ISM (e. g. ISlyz et al1l2005l: 
'Joung & M ac Lowl l2006l: 'Joung et al." 2009), the spatial reso- 
lutions are usually above 1 pc, and in simu lations of galactic 
disks, the resolutio n is even coarser (e. g. lAgertz et^ I2QQ9I: 
iTasker & Tai]l2009h . Thus, as shown here, it is quite likely that 
the mass fractions of the gas in the diff'erent phases would change 
significantly if smaller length scales could be resolved, particu- 
larly in the cold phase and in the low-temperature segment of 
the unstable phase. This has a potential impact on calculations 
of the star formation rate, because only cold gas can turn into 
stars. A subgrid-model for star formation that takes the multi- 
phase structure of the I SM into account has already been de- 
veloped and applied by ISpringel & Hernquisi (12003.) , although 
these authors omitted the influence of turbulent motions that we 
focussed on. However, our results indicate that the errors due 
to unresolved thermodynamic processes will be small if the gas 
density in the grid cells where sink particles representing star- 
forming clouds are added is much higher than 1.0 cm"^. Only in 
regions that just begin to contract, a correction of the star forma- 
tion efficiency by the mass fraction of cold gas might slow down 
the growth of the stellar mass in disk-galaxy simulations. In the 
dense clumps that are formed subsequently, however, most of the 
gas would cool down to molecular-cloud temperatures if a suflft- 



ciently low temperature floor was applied, and it can be assumed 
that the local star formation efficie ncy is regulated mai nly by 
gravity and small-scale turbulence ((&um holz & McKe^ [20051: 
iPadoan & Nordlundll2009 l : iKlessen & HennebeUe 20 10). 

For further studies of turbulence in the multi-phase ISM, 
adaptive mesh refinement (AMR) is one method to in- 
crea se the numerica l resolution in such simulations. Although 
Nikl aus et a n (I2009h show that there are serious difl&culties in 
reproducing the properties of 2-dimensional colliding flows by 
means of AMR even if sophisticated refinement criteria are ap- 
plied, we expect that the situation will be more favourable to 
3 -dimensional simulations of forced turbulence. If the low en- 
strophy (i. e., the volume-integrated square of the vorticity) of 
compact cold gas regions could be confirmed in AMR sim- 
ulations, it would point to an additional driving mechanism 
to sustain supersonic turbulence in molecular cloud interiors. 
Examples include external driving by gravitational accretion 
(e. g. Klessen & Hennebelle 2010) or internal driving by mas- 
sive stars (e. g. iFall et al .1120101) . Moreover, there is also room 
for methodological improvements in the cooling. Firstly, the 
range of the cooling function should be extended to temperatures 
higher than 10 000 K. On the other hand, the cooling of the gas at 
low temperatures will be altered if chemical processes, in partic- 
ular, the generation and destruction of molecular hydrogen, are 
explicitly included (iGlover et al.ll2QlQl) . Secondly, it is not clear 
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whether an explicit scheme with subcycling will be sufficient to 
treat the source terms in the energy equation if the effective nu- 
merical resolution becomes very high. For subsequent simula- 
tions it would also be interesting to study the effect of magnetic 
fields. 
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